module mod_clm_cnallocation
#ifdef CN
  !
  ! Module holding routines used in allocation model for coupled carbon
  ! nitrogen code.
  !
  use mod_intkinds
  use mod_realkinds
  use mod_runparams, only : dtsrf
  use mod_clm_varcon, only : dzsoi_decomp
  use mod_stdio
  use mod_mpmessage
  use mod_clm_varctl, only : use_c13, use_c14

  implicit none

  private

  save

  public :: CNAllocationInit         ! Initialization
  public :: CNAllocation             ! run method
  public :: CNAllocation_Carbon_only ! Return Carbon_only status

  ! Supplemental Nitrogen for all PFT's
  character(len=*), parameter, public :: suplnAll = 'ALL'
  ! No supplemental Nitrogen
  character(len=*), parameter, public :: suplnNon = 'NONE'
  ! Supplemental Nitrogen mode
  character(len=15), public :: suplnitro = suplnNon
  ! Carbon only mode
  ! (Nitrogen is prescribed NOT prognostic)
  logical, public :: carbon_only = .false.

  real(rk8) :: dt           !decomp timestep (seconds)
  real(rk8) :: bdnr         !bulk denitrification rate (1/s)
  real(rk8) :: dayscrecover !number of days to recover negative cpool
  !reproduction allocation coefficient
  real(rk8), pointer, contiguous, dimension(:) :: arepr
  !root allocation coefficient
  real(rk8), pointer, contiguous, dimension(:) :: aroot
  !column-level plant N demand
  real(rk8), pointer, contiguous, dimension(:) :: col_plant_ndemand
  ! Prognostic crop receives supplemental Nitrogen
  ! logical :: crop_supln  = .false.

  contains

  logical function CNAllocation_Carbon_only()
    CNAllocation_Carbon_only = carbon_only
  end function CNAllocation_Carbon_only

  subroutine CNAllocationInit ( lbc, ubc, lbp, ubp )
    use mod_clm_varcon, only : secspday
    use mod_clm_surfrd, only : crop_prog
    implicit none
    integer(ik4), intent(in) :: lbc, ubc        ! column-index bounds
    integer(ik4), intent(in) :: lbp, ubp        ! pft-index bounds
    character(len=32) :: subname = 'CNAllocationInit'
    if ( crop_prog )then
      allocate(arepr(lbp:ubp))
      allocate(aroot(lbp:ubp))
      arepr(:) = nan
      aroot(:) = nan
    end if
    allocate(col_plant_ndemand(lbc:ubc))
    col_plant_ndemand(:) = nan

    ! set time steps
    dt = dtsrf

    ! set some space-and-time constant parameters
    bdnr         = 0.5_rk8 * (dt/secspday)
    dayscrecover = 30.0_rk8

    ! Change namelist settings into private logical variables
    select case(suplnitro)
      case(suplnNon)
        carbon_only = .false.
      case(suplnAll)
        carbon_only = .true.
      case default
        write(stderr,*) &
                'Supplemental Nitrogen flag (suplnitro) can only be: ', &
               suplnNon, ' or ', suplnAll
        call fatal(__FILE__,__LINE__, &
          trim(subname)//'ERROR: supplemental Nitrogen flag is not correct' )
    end select
  end subroutine CNAllocationInit

  subroutine CNAllocation (lbp, ubp, lbc, ubc, &
       num_soilc, filter_soilc, num_soilp, filter_soilp )
    use mod_clm_type
    use mod_clm_subgridave, only : p2c
    use mod_clm_varpar, only : nlevsoi, nlevdecomp
#ifdef NITRIF_DENITRIF
    use mod_clm_varcon, only : nitrif_n2o_loss_frac
#endif
    use mod_clm_pftvarcon, only : npcropmin, declfact, bfact, aleaff,  &
            arootf, astemf, arooti, fleafi, allconsl, allconss, grperc, &
            grpnow, nsoybean
    use mod_clm_varcon, only : secspday, istsoil, istcrop
    use mod_clm_varpar, only : max_pft_per_col
    !use mod_clm_pftvarcon, only : nbrdlf_evr_trp_tree

    implicit none
    integer(ik4), intent(in) :: lbp, ubp  ! pft-index bounds
    integer(ik4), intent(in) :: lbc, ubc  ! column-index bounds
    ! number of soil columns in filter
    integer(ik4), intent(in) :: num_soilc
    ! filter for soil columns
    integer(ik4), intent(in) :: filter_soilc(ubc-lbc+1)
    integer(ik4), intent(in) :: num_soilp       ! number of soil pfts in filter
    integer(ik4), intent(in) :: filter_soilp(ubp-lbp+1) ! filter for soil pfts

    ! pft level
    integer(ik4), pointer, contiguous :: ivt(:)        ! pft vegetation type
    integer(ik4), pointer, contiguous :: pcolumn(:)    ! pft's column index
    integer(ik4), pointer, contiguous :: pfti(:)       ! initial pft index in landunit
    real(rk8), pointer, contiguous :: lgsf(:)       ! long growing season factor [0-1]
    real(rk8), pointer, contiguous :: xsmrpool(:)   ! (gC/m2) temporary photosynthate C pool
    real(rk8), pointer, contiguous :: retransn(:)   ! (gN/m2) plant pool of retranslocated N
    ! sunlit leaf-level photosynthesis (umol CO2 /m**2/ s)
    real(rk8), pointer, contiguous :: psnsun(:)
    ! shaded leaf-level photosynthesis (umol CO2 /m**2/ s)
    real(rk8), pointer, contiguous :: psnsha(:)

    ! C13 sunlit leaf-level photosynthesis (umol CO2 /m**2/ s)
    real(rk8), pointer, contiguous :: c13_psnsun(:)
    ! C13 shaded leaf-level photosynthesis (umol CO2 /m**2/ s)
    real(rk8), pointer, contiguous :: c13_psnsha(:)

    ! C14 sunlit leaf-level photosynthesis (umol CO2 /m**2/ s)
    real(rk8), pointer, contiguous :: c14_psnsun(:)
    ! C14 shaded leaf-level photosynthesis (umol CO2 /m**2/ s)
    real(rk8), pointer, contiguous :: c14_psnsha(:)

    real(rk8), pointer, contiguous :: laisun(:)     ! sunlit projected leaf area index
    real(rk8), pointer, contiguous :: laisha(:)     ! shaded projected leaf area index
    real(rk8), pointer, contiguous :: leafc(:)
    real(rk8), pointer, contiguous :: frootc(:)
    real(rk8), pointer, contiguous :: livestemc(:)
    real(rk8), pointer, contiguous :: leaf_mr(:)
    real(rk8), pointer, contiguous :: froot_mr(:)
    real(rk8), pointer, contiguous :: livestem_mr(:)
    real(rk8), pointer, contiguous :: livecroot_mr(:)
    real(rk8), pointer, contiguous :: grain_mr(:)
    real(rk8), pointer, contiguous :: leaf_curmr(:)
    real(rk8), pointer, contiguous :: froot_curmr(:)
    real(rk8), pointer, contiguous :: livestem_curmr(:)
    real(rk8), pointer, contiguous :: livecroot_curmr(:)
    real(rk8), pointer, contiguous :: grain_curmr(:)
    real(rk8), pointer, contiguous :: leaf_xsmr(:)
    real(rk8), pointer, contiguous :: froot_xsmr(:)
    real(rk8), pointer, contiguous :: livestem_xsmr(:)
    real(rk8), pointer, contiguous :: livecroot_xsmr(:)
    real(rk8), pointer, contiguous :: grain_xsmr(:)
    ! column level
    real(rk8), pointer, contiguous :: sminn_vr(:,:)   ! (gN/m3) soil mineral N
    ! ecophysiological constants
    ! binary flag for woody lifeform (1=woody, 0=not woody)
    real(rk8), pointer, contiguous :: woody(:)
    ! allocation parameter: new fine root C per new leaf C (gC/gC)
    real(rk8), pointer, contiguous :: froot_leaf(:)
    ! allocation parameter: new coarse root C per new stem C (gC/gC)
    real(rk8), pointer, contiguous :: croot_stem(:)
    ! allocation parameter: new stem c per new leaf C (gC/gC)
    real(rk8), pointer, contiguous :: stem_leaf(:)
    ! allocation parameter: fraction of new wood that is live
    ! (phloem and ray parenchyma) (no units)
    real(rk8), pointer, contiguous :: flivewd(:)
    real(rk8), pointer, contiguous :: leafcn(:)   ! leaf C:N (gC/gN)
    real(rk8), pointer, contiguous :: frootcn(:)  ! fine root C:N (gC/gN)
    ! live wood (phloem and ray parenchyma) C:N (gC/gN)
    real(rk8), pointer, contiguous :: livewdcn(:)
    ! dead wood (xylem and heartwood) C:N (gC/gN)
    real(rk8), pointer, contiguous :: deadwdcn(:)
    ! allocation parameter: fraction of allocation that goes to currently
    ! displayed growth, remainder to storage
    real(rk8), pointer, contiguous :: fcur2(:)
    integer(ik4), pointer, contiguous :: plandunit(:) ! index into landunit level quantities
    integer(ik4), pointer, contiguous :: clandunit(:) ! index into landunit level quantities
    integer(ik4), pointer, contiguous :: itypelun(:) ! landunit type
    logical, pointer, contiguous :: croplive(:)     ! flag, true if planted, not harvested
    integer(ik4), pointer, contiguous :: peaklai(:) ! 1: max allowed lai; 0: not at max
    real(rk8), pointer, contiguous :: gddmaturity(:) ! gdd needed to harvest
    ! heat unit index needed from planting to leaf emergence
    real(rk8), pointer, contiguous :: huileaf(:)
    real(rk8), pointer, contiguous :: huigrain(:) ! same to reach vegetative maturity
    real(rk8), pointer, contiguous :: hui(:)      ! =gdd since planting (gddplant)
    real(rk8), pointer, contiguous :: leafout(:)  ! =gdd from top soil layer temperature
    real(rk8), pointer, contiguous :: aleafi(:) ! saved allocation coefficient from phase 2
    real(rk8), pointer, contiguous :: astemi(:) ! saved allocation coefficient from phase 2
    real(rk8), pointer, contiguous :: aleaf(:)      ! leaf allocation coefficient
    real(rk8), pointer, contiguous :: astem(:)      ! stem allocation coefficient
    real(rk8), pointer, contiguous :: graincn(:)    ! grain C:N (gC/gN)
    real(rk8), pointer, contiguous :: fleafcn(:)    ! leaf c:n during organ fill
    real(rk8), pointer, contiguous :: fstemcn(:)    ! stem c:n during organ fill
    real(rk8), pointer, contiguous :: ffrootcn(:)   ! froot c:n during organ fill

    ! pft level
    real(rk8), pointer, contiguous :: grain_flag(:) ! 1: grain fill stage; 0: not
    real(rk8), pointer, contiguous :: gpp(:)     ! GPP flux before downregulation (gC/m2/s)
    real(rk8), pointer, contiguous :: availc(:)  ! C flux available for allocation (gC/m2/s)
    ! C flux assigned to recovery of negative cpool (gC/m2/s)
    real(rk8), pointer, contiguous :: xsmrpool_recover(:)
    real(rk8), pointer, contiguous :: c_allometry(:)        ! C allocation index (DIM)
    real(rk8), pointer, contiguous :: n_allometry(:)        ! N allocation index (DIM)
    ! N flux required to support initial GPP (gN/m2/s)
    real(rk8), pointer, contiguous :: plant_ndemand(:)
    ! temporary annual sum of potential GPP
    real(rk8), pointer, contiguous :: tempsum_potential_gpp(:)
    ! temporary annual max of retranslocated N pool (gN/m2)
    real(rk8), pointer, contiguous :: tempmax_retransn(:)
    ! annual sum of potential GPP
    real(rk8), pointer, contiguous :: annsum_potential_gpp(:)
    ! N flux available from retranslocation pool (gN/m2/s)
    real(rk8), pointer, contiguous :: avail_retransn(:)
    ! annual max of retranslocated N pool
    real(rk8), pointer, contiguous :: annmax_retransn(:)
    ! total allocated N flux (gN/m2/s)
    real(rk8), pointer, contiguous :: plant_nalloc(:)
    ! total allocated C flux (gC/m2/s)
    real(rk8), pointer, contiguous :: plant_calloc(:)
    ! C flux not allocated due to downregulation (gC/m2/s)
    real(rk8), pointer, contiguous :: excess_cflux(:)
    ! fractional reduction in GPP due to N limitation (DIM)
    real(rk8), pointer, contiguous :: downreg(:)
    ! annual sum of NPP, for wood allocation
    real(rk8), pointer, contiguous :: annsum_npp(:)
    real(rk8), pointer, contiguous :: cpool_to_xsmrpool(:)
    real(rk8), pointer, contiguous :: psnsun_to_cpool(:)
    real(rk8), pointer, contiguous :: psnshade_to_cpool(:)

    real(rk8), pointer, contiguous :: c13_psnsun_to_cpool(:)
    real(rk8), pointer, contiguous :: c13_psnshade_to_cpool(:)

    real(rk8), pointer, contiguous :: c14_psnsun_to_cpool(:)
    real(rk8), pointer, contiguous :: c14_psnshade_to_cpool(:)

    real(rk8), pointer, contiguous :: cpool_to_leafc(:)
    real(rk8), pointer, contiguous :: cpool_to_leafc_storage(:)
    real(rk8), pointer, contiguous :: cpool_to_frootc(:)
    real(rk8), pointer, contiguous :: cpool_to_frootc_storage(:)
    real(rk8), pointer, contiguous :: cpool_to_livestemc(:)
    real(rk8), pointer, contiguous :: cpool_to_livestemc_storage(:)
    real(rk8), pointer, contiguous :: cpool_to_deadstemc(:)
    real(rk8), pointer, contiguous :: cpool_to_deadstemc_storage(:)
    real(rk8), pointer, contiguous :: cpool_to_livecrootc(:)
    real(rk8), pointer, contiguous :: cpool_to_livecrootc_storage(:)
    real(rk8), pointer, contiguous :: cpool_to_deadcrootc(:)
    real(rk8), pointer, contiguous :: cpool_to_deadcrootc_storage(:)
    ! allocation to growth respiration storage (gC/m2/s)
    real(rk8), pointer, contiguous :: cpool_to_gresp_storage(:)
    ! deployment of retranslocated N (gN/m2/s)
    real(rk8), pointer, contiguous :: retransn_to_npool(:)
    ! deployment of soil mineral N uptake (gN/m2/s)
    real(rk8), pointer, contiguous :: sminn_to_npool(:)
    ! allocation to grain C (gC/m2/s)
    real(rk8), pointer, contiguous :: cpool_to_grainc(:)
    ! allocation to grain C storage (gC/m2/s)
    real(rk8), pointer, contiguous :: cpool_to_grainc_storage(:)
    ! allocation to grain N (gN/m2/s)
    real(rk8), pointer, contiguous :: npool_to_grainn(:)
    ! allocation to grain N storage (gN/m2/s)
    real(rk8), pointer, contiguous :: npool_to_grainn_storage(:)
    ! allocation to leaf N (gN/m2/s)
    real(rk8), pointer, contiguous :: npool_to_leafn(:)
    ! allocation to leaf N storage (gN/m2/s)
    real(rk8), pointer, contiguous :: npool_to_leafn_storage(:)
    ! allocation to fine root N (gN/m2/s)
    real(rk8), pointer, contiguous :: npool_to_frootn(:)
    ! allocation to fine root N storage (gN/m2/s)
    real(rk8), pointer, contiguous :: npool_to_frootn_storage(:)
    real(rk8), pointer, contiguous :: npool_to_livestemn(:)
    real(rk8), pointer, contiguous :: npool_to_livestemn_storage(:)
    real(rk8), pointer, contiguous :: npool_to_deadstemn(:)
    real(rk8), pointer, contiguous :: npool_to_deadstemn_storage(:)
    real(rk8), pointer, contiguous :: npool_to_livecrootn(:)
    real(rk8), pointer, contiguous :: npool_to_livecrootn_storage(:)
    real(rk8), pointer, contiguous :: npool_to_deadcrootn(:)
    real(rk8), pointer, contiguous :: npool_to_deadcrootn_storage(:)
    ! column level
    ! fraction of potential immobilization (no units)
    real(rk8), pointer, contiguous :: fpi(:)
    ! fraction of potential gpp (no units)
    real(rk8), pointer, contiguous :: fpg(:)
    real(rk8), pointer, contiguous :: potential_immob(:)
    real(rk8), pointer, contiguous :: actual_immob(:)
    real(rk8), pointer, contiguous :: sminn_to_plant(:)
    ! fraction of potential immobilization (no units)
    real(rk8), pointer, contiguous :: fpi_vr(:,:)
    real(rk8), pointer, contiguous, dimension(:,:) :: h2osoi_vol
#ifndef NITRIF_DENITRIF
    real(rk8), pointer, contiguous :: sminn_to_denit_excess_vr(:,:)
#else
    ! (unitless) relative compettiveness of plants for NO3
    real(rk8), parameter :: compet_plant_no3   = 1.0
    ! (unitless) relative compettiveness of plants for NH4
    real(rk8), parameter :: compet_plant_nh4   = 1.0
    ! (unitless) relative competitiveness of immobilizers for NO3
    real(rk8), parameter :: compet_decomp_no3  = 1.0
    ! (unitless) relative competitiveness of immobilizers for NH4
    real(rk8), parameter :: compet_decomp_nh4  = 1.0
    ! (unitless) relative competitiveness of denitrifiers for NO3
    real(rk8), parameter :: compet_denit       = 1.0
    ! (unitless) relative competitiveness of nitrifiers for NH4
    real(rk8), parameter :: compet_nit         = 1.0
    ! fraction of potential immobilization supplied by no3(no units)
    real(rk8) :: fpi_no3_vr(lbc:ubc,1:nlevdecomp)
    ! fraction of potential immobilization supplied by nh4 (no units)
    real(rk8) :: fpi_nh4_vr(lbc:ubc,1:nlevdecomp)
    real(rk8) :: sum_nh4_demand(lbc:ubc,1:nlevdecomp), &
                 sum_nh4_demand_scaled(lbc:ubc,1:nlevdecomp)
    real(rk8) :: sum_no3_demand(lbc:ubc,1:nlevdecomp), &
                 sum_no3_demand_scaled(lbc:ubc,1:nlevdecomp)
    real(rk8), pointer, contiguous :: smin_no3_vr(:,:) ! (gN/m3) soil mineral NO3
    real(rk8), pointer, contiguous :: smin_nh4_vr(:,:) ! (gN/m3) soil mineral NH4
    real(rk8), pointer, contiguous :: f_nit_vr(:,:)    ! (gN/m3/s) soil nitrification flux
    real(rk8), pointer, contiguous :: f_denit_vr(:,:)  ! (gN/m3/s) soil denitrification flux
    ! (gN/m3/s) potential soil nitrification flux
    real(rk8), pointer, contiguous :: pot_f_nit_vr(:,:)
    ! (gN/m3/s) potential soil denitrification flux
    real(rk8), pointer, contiguous :: pot_f_denit_vr(:,:)
    real(rk8), pointer, contiguous :: actual_immob_no3_vr(:,:)
    real(rk8), pointer, contiguous :: actual_immob_nh4_vr(:,:)
    real(rk8), pointer, contiguous :: smin_no3_to_plant_vr(:,:)
    real(rk8), pointer, contiguous :: smin_nh4_to_plant_vr(:,:)
    ! ratio of N2 to N2O production by denitrification [gN/gN]
    real(rk8), pointer, contiguous :: n2_n2o_ratio_denit_vr(:,:)
    ! flux of N2O from denitrification [gN/m3/s]
    real(rk8), pointer, contiguous :: f_n2o_denit_vr(:,:)
    ! flux of N2O from nitrification [gN/m3/s]
    real(rk8), pointer, contiguous :: f_n2o_nit_vr(:,:)
    real(rk8), pointer, contiguous :: f_n2o_tot_vr(:,:)
    real(rk8), pointer, contiguous :: fr_WFPS(:,:)
#endif
    real(rk8), pointer, contiguous :: sminn_to_plant_vr(:,:)
    real(rk8), pointer, contiguous :: supplement_to_sminn_vr(:,:)
    real(rk8), pointer, contiguous :: nfixation_prof(:,:)
    real(rk8), pointer, contiguous :: potential_immob_vr(:,:)
    real(rk8), pointer, contiguous :: actual_immob_vr(:,:)

    real(rk8), pointer, contiguous :: leafn_to_retransn(:)
    real(rk8), pointer, contiguous :: frootn_to_retransn(:)
    real(rk8), pointer, contiguous :: livestemn_to_retransn(:)

    integer(ik4) :: c,p,l  !indices
    integer(ik4) :: fp     !lake filter pft index
    integer(ik4) :: fc     !lake filter column index
    real(rk8):: mr         !maintenance respiration (gC/m2/s)
    real(rk8):: f1,f2,f3,f4,g1,g2  !allocation parameters
    real(rk8):: cnl,cnfr,cnlw,cndw !C:N ratios for leaf, fine root, and wood
    real(rk8):: fcur           !fraction of current psn displayed as growth
    real(rk8):: gresp_storage  !temporary variable for growth resp to storage
    !temporary variable for total new leaf carbon allocation
    real(rk8):: nlc
    real(rk8):: curmr, curmr_ratio   !xsmrpool temporary variables
    !total column N demand (gN/m3/s) at a given level
    real(rk8):: sum_ndemand_vr(lbc:ubc, 1:nlevdecomp)
    integer(ik4) :: j     ! level index
    real(rk8):: nuptake_prof(lbc:ubc, 1:nlevdecomp)
    real(rk8):: sminn_tot(lbc:ubc)
    real(rk8) f5          !grain allocation parameter
    real(rk8) cng         !C:N ratio for grain (= cnlw for now; slevis)
    real(rk8) fleaf       !fraction allocated to leaf
    real(rk8) t1          !temporary variable

#ifndef NITRIF_DENITRIF
    !flag for N limitation
    integer(ik4) :: nlimit(lbc:ubc,0:nlevdecomp)
    real(rk8):: residual_sminn_vr(lbc:ubc, 1:nlevdecomp)
    real(rk8):: residual_sminn(lbc:ubc)
#else
    !flag for NO3 limitation
    integer(ik4) :: nlimit_no3(lbc:ubc,0:nlevdecomp)
    !flag for NH4 limitation
    integer(ik4) :: nlimit_nh4(lbc:ubc,0:nlevdecomp)
    real(rk8):: residual_smin_nh4_vr(lbc:ubc, 1:nlevdecomp)
    real(rk8):: residual_smin_no3_vr(lbc:ubc, 1:nlevdecomp)
    real(rk8):: residual_smin_nh4(lbc:ubc)
    real(rk8):: residual_smin_no3(lbc:ubc)
#endif
    real(rk8):: residual_plant_ndemand(lbc:ubc)
    logical :: iswood(lbp:ubp)

    ! Assign local pointers to derived type arrays (in)
    ivt                         => clm3%g%l%c%p%itype
    pcolumn                     => clm3%g%l%c%p%column
    plandunit                   => clm3%g%l%c%p%landunit
    clandunit                   => clm3%g%l%c%landunit
    pfti                        => clm3%g%l%c%pfti
    itypelun                    => clm3%g%l%itype
    lgsf                        => clm3%g%l%c%p%pepv%lgsf
    xsmrpool                    => clm3%g%l%c%p%pcs%xsmrpool
    retransn                    => clm3%g%l%c%p%pns%retransn
    psnsun                      => clm3%g%l%c%p%pcf%psnsun
    psnsha                      => clm3%g%l%c%p%pcf%psnsha
    h2osoi_vol                  => clm3%g%l%c%cws%h2osoi_vol
    if ( use_c13 ) then
      c13_psnsun                  => clm3%g%l%c%p%pc13f%psnsun
      c13_psnsha                  => clm3%g%l%c%p%pc13f%psnsha
      c13_psnsun_to_cpool         => clm3%g%l%c%p%pc13f%psnsun_to_cpool
      c13_psnshade_to_cpool       => clm3%g%l%c%p%pc13f%psnshade_to_cpool
    end if
    if ( use_c14 ) then
      c14_psnsun                  => clm3%g%l%c%p%pc14f%psnsun
      c14_psnsha                  => clm3%g%l%c%p%pc14f%psnsha
      c14_psnsun_to_cpool         => clm3%g%l%c%p%pc14f%psnsun_to_cpool
      c14_psnshade_to_cpool       => clm3%g%l%c%p%pc14f%psnshade_to_cpool
    end if
    laisun                      => clm3%g%l%c%p%pps%laisun
    laisha                      => clm3%g%l%c%p%pps%laisha
    leafc                       => clm3%g%l%c%p%pcs%leafc
    frootc                      => clm3%g%l%c%p%pcs%frootc
    livestemc                   => clm3%g%l%c%p%pcs%livestemc
    leaf_mr                     => clm3%g%l%c%p%pcf%leaf_mr
    froot_mr                    => clm3%g%l%c%p%pcf%froot_mr
    livestem_mr                 => clm3%g%l%c%p%pcf%livestem_mr
    livecroot_mr                => clm3%g%l%c%p%pcf%livecroot_mr
    grain_mr                    => clm3%g%l%c%p%pcf%grain_mr
    leaf_curmr                  => clm3%g%l%c%p%pcf%leaf_curmr
    froot_curmr                 => clm3%g%l%c%p%pcf%froot_curmr
    livestem_curmr              => clm3%g%l%c%p%pcf%livestem_curmr
    livecroot_curmr             => clm3%g%l%c%p%pcf%livecroot_curmr
    grain_curmr                 => clm3%g%l%c%p%pcf%grain_curmr
    leaf_xsmr                   => clm3%g%l%c%p%pcf%leaf_xsmr
    froot_xsmr                  => clm3%g%l%c%p%pcf%froot_xsmr
    livestem_xsmr               => clm3%g%l%c%p%pcf%livestem_xsmr
    livecroot_xsmr              => clm3%g%l%c%p%pcf%livecroot_xsmr
    grain_xsmr                  => clm3%g%l%c%p%pcf%grain_xsmr
    sminn_vr                       => clm3%g%l%c%cns%sminn_vr
    woody                       => pftcon%woody
    froot_leaf                  => pftcon%froot_leaf
    croot_stem                  => pftcon%croot_stem
    stem_leaf                   => pftcon%stem_leaf
    flivewd                     => pftcon%flivewd
    leafcn                      => pftcon%leafcn
    frootcn                     => pftcon%frootcn
    livewdcn                    => pftcon%livewdcn
    deadwdcn                    => pftcon%deadwdcn
    fcur2                       => pftcon%fcur
    gddmaturity                 => clm3%g%l%c%p%pps%gddmaturity
    huileaf                     => clm3%g%l%c%p%pps%huileaf
    huigrain                    => clm3%g%l%c%p%pps%huigrain
    hui                         => clm3%g%l%c%p%pps%gddplant
    leafout                     => clm3%g%l%c%p%pps%gddtsoi
    croplive                    => clm3%g%l%c%p%pps%croplive
    peaklai                     => clm3%g%l%c%p%pps%peaklai
    graincn                     => pftcon%graincn
    fleafcn                     => pftcon%fleafcn
    ffrootcn                    => pftcon%ffrootcn
    fstemcn                     => pftcon%fstemcn
    ! Assign local pointers to derived type arrays (out)
    grain_flag                  => clm3%g%l%c%p%pepv%grain_flag
    gpp                         => clm3%g%l%c%p%pepv%gpp
    availc                      => clm3%g%l%c%p%pepv%availc
    xsmrpool_recover            => clm3%g%l%c%p%pepv%xsmrpool_recover
    c_allometry                 => clm3%g%l%c%p%pepv%c_allometry
    n_allometry                 => clm3%g%l%c%p%pepv%n_allometry
    plant_ndemand               => clm3%g%l%c%p%pepv%plant_ndemand
    tempsum_potential_gpp       => clm3%g%l%c%p%pepv%tempsum_potential_gpp
    tempmax_retransn            => clm3%g%l%c%p%pepv%tempmax_retransn
    annsum_potential_gpp        => clm3%g%l%c%p%pepv%annsum_potential_gpp
    avail_retransn              => clm3%g%l%c%p%pepv%avail_retransn
    annmax_retransn             => clm3%g%l%c%p%pepv%annmax_retransn
    plant_nalloc                => clm3%g%l%c%p%pepv%plant_nalloc
    plant_calloc                => clm3%g%l%c%p%pepv%plant_calloc
    excess_cflux                => clm3%g%l%c%p%pepv%excess_cflux
    downreg                     => clm3%g%l%c%p%pepv%downreg
    annsum_npp                  => clm3%g%l%c%p%pepv%annsum_npp
    cpool_to_xsmrpool           => clm3%g%l%c%p%pcf%cpool_to_xsmrpool
    psnsun_to_cpool             => clm3%g%l%c%p%pcf%psnsun_to_cpool
    psnshade_to_cpool           => clm3%g%l%c%p%pcf%psnshade_to_cpool
    cpool_to_leafc              => clm3%g%l%c%p%pcf%cpool_to_leafc
    cpool_to_leafc_storage      => clm3%g%l%c%p%pcf%cpool_to_leafc_storage
    cpool_to_frootc             => clm3%g%l%c%p%pcf%cpool_to_frootc
    cpool_to_frootc_storage     => clm3%g%l%c%p%pcf%cpool_to_frootc_storage
    cpool_to_livestemc          => clm3%g%l%c%p%pcf%cpool_to_livestemc
    cpool_to_livestemc_storage  => clm3%g%l%c%p%pcf%cpool_to_livestemc_storage
    cpool_to_deadstemc          => clm3%g%l%c%p%pcf%cpool_to_deadstemc
    cpool_to_deadstemc_storage  => clm3%g%l%c%p%pcf%cpool_to_deadstemc_storage
    cpool_to_livecrootc         => clm3%g%l%c%p%pcf%cpool_to_livecrootc
    cpool_to_livecrootc_storage => clm3%g%l%c%p%pcf%cpool_to_livecrootc_storage
    cpool_to_deadcrootc         => clm3%g%l%c%p%pcf%cpool_to_deadcrootc
    cpool_to_deadcrootc_storage => clm3%g%l%c%p%pcf%cpool_to_deadcrootc_storage
    cpool_to_gresp_storage      => clm3%g%l%c%p%pcf%cpool_to_gresp_storage
    cpool_to_grainc             => clm3%g%l%c%p%pcf%cpool_to_grainc
    cpool_to_grainc_storage     => clm3%g%l%c%p%pcf%cpool_to_grainc_storage
    npool_to_grainn             => clm3%g%l%c%p%pnf%npool_to_grainn
    npool_to_grainn_storage     => clm3%g%l%c%p%pnf%npool_to_grainn_storage
    retransn_to_npool           => clm3%g%l%c%p%pnf%retransn_to_npool
    sminn_to_npool              => clm3%g%l%c%p%pnf%sminn_to_npool
    npool_to_leafn              => clm3%g%l%c%p%pnf%npool_to_leafn
    npool_to_leafn_storage      => clm3%g%l%c%p%pnf%npool_to_leafn_storage
    npool_to_frootn             => clm3%g%l%c%p%pnf%npool_to_frootn
    npool_to_frootn_storage     => clm3%g%l%c%p%pnf%npool_to_frootn_storage
    npool_to_livestemn          => clm3%g%l%c%p%pnf%npool_to_livestemn
    npool_to_livestemn_storage  => clm3%g%l%c%p%pnf%npool_to_livestemn_storage
    npool_to_deadstemn          => clm3%g%l%c%p%pnf%npool_to_deadstemn
    npool_to_deadstemn_storage  => clm3%g%l%c%p%pnf%npool_to_deadstemn_storage
    npool_to_livecrootn         => clm3%g%l%c%p%pnf%npool_to_livecrootn
    npool_to_livecrootn_storage => clm3%g%l%c%p%pnf%npool_to_livecrootn_storage
    npool_to_deadcrootn         => clm3%g%l%c%p%pnf%npool_to_deadcrootn
    npool_to_deadcrootn_storage => clm3%g%l%c%p%pnf%npool_to_deadcrootn_storage
    leafn_to_retransn           => clm3%g%l%c%p%pnf%leafn_to_retransn
    frootn_to_retransn          => clm3%g%l%c%p%pnf%frootn_to_retransn
    livestemn_to_retransn       => clm3%g%l%c%p%pnf%livestemn_to_retransn
    fpg                         => clm3%g%l%c%cps%fpg
    potential_immob             => clm3%g%l%c%cnf%potential_immob
    actual_immob                => clm3%g%l%c%cnf%actual_immob
    sminn_to_plant              => clm3%g%l%c%cnf%sminn_to_plant
    fpi                     => clm3%g%l%c%cps%fpi
    fpi_vr                  => clm3%g%l%c%cps%fpi_vr
#ifndef NITRIF_DENITRIF
    sminn_to_denit_excess_vr => clm3%g%l%c%cnf%sminn_to_denit_excess_vr
#else
    pot_f_nit_vr           => clm3%g%l%c%cnf%pot_f_nit_vr
    pot_f_denit_vr         => clm3%g%l%c%cnf%pot_f_denit_vr
    f_nit_vr               => clm3%g%l%c%cnf%f_nit_vr
    f_denit_vr             => clm3%g%l%c%cnf%f_denit_vr
    actual_immob_no3_vr    => clm3%g%l%c%cnf%actual_immob_no3_vr
    actual_immob_nh4_vr    => clm3%g%l%c%cnf%actual_immob_nh4_vr
    smin_no3_to_plant_vr   => clm3%g%l%c%cnf%smin_no3_to_plant_vr
    smin_nh4_to_plant_vr   => clm3%g%l%c%cnf%smin_nh4_to_plant_vr
    n2_n2o_ratio_denit_vr  => clm3%g%l%c%cnf%n2_n2o_ratio_denit_vr
    f_n2o_denit_vr         => clm3%g%l%c%cnf%f_n2o_denit_vr
    f_n2o_nit_vr           => clm3%g%l%c%cnf%f_n2o_nit_vr
    f_n2o_tot_vr           => clm3%g%l%c%cnf%f_n2o_tot_vr
    fr_WFPS                => clm3%g%l%c%cnf%fr_WFPS
    smin_no3_vr            => clm3%g%l%c%cns%smin_no3_vr
    smin_nh4_vr            => clm3%g%l%c%cns%smin_nh4_vr
#endif
    supplement_to_sminn_vr => clm3%g%l%c%cnf%supplement_to_sminn_vr
    sminn_to_plant_vr      => clm3%g%l%c%cnf%sminn_to_plant_vr
    nfixation_prof         => clm3%g%l%c%cps%nfixation_prof
    potential_immob_vr     => clm3%g%l%c%cnf%potential_immob_vr
    actual_immob_vr        => clm3%g%l%c%cnf%actual_immob_vr
    aleafi                 => clm3%g%l%c%p%pps%aleafi
    astemi                 => clm3%g%l%c%p%pps%astemi
    aleaf                  => clm3%g%l%c%p%pps%aleaf
    astem                  => clm3%g%l%c%p%pps%astem

    ! set time steps
    dt = dtsrf

    do p = lbp, ubp
      iswood(p) = abs(woody(ivt(p))-1._rk8) < epsilon(1.0)
    end do

    ! loop over pfts to assess the total plant N demand
    do fp = 1, num_soilp
      p = filter_soilp(fp)

      ! get the time step total gross photosynthesis
      ! this is coming from the canopy fluxes code, and is the
      ! gpp that is used to control stomatal conductance.
      ! For the nitrogen downregulation code, this is assumed
      ! to be the potential gpp, and the actual gpp will be
      ! reduced due to N limitation.

      ! Convert psn from umol/m2/s -> gC/m2/s

      ! The input psn (psnsun and psnsha) are expressed per unit LAI
      ! in the sunlit and shaded canopy, respectively. These need to be
      ! scaled by laisun and laisha to get the total gpp for allocation

      psnsun_to_cpool(p) = psnsun(p) * laisun(p) * 12.011e-6_rk8
      psnshade_to_cpool(p) = psnsha(p) * laisha(p) * 12.011e-6_rk8

      if ( use_c13 ) then
        c13_psnsun_to_cpool(p) = c13_psnsun(p) * laisun(p) * 12.011e-6_rk8
        c13_psnshade_to_cpool(p) = c13_psnsha(p) * laisha(p) * 12.011e-6_rk8
      end if

      if ( use_c14 ) then
        c14_psnsun_to_cpool(p) = c14_psnsun(p) * laisun(p) * 12.011e-6_rk8
        c14_psnshade_to_cpool(p) = c14_psnsha(p) * laisha(p) * 12.011e-6_rk8
      end if

      gpp(p) = psnsun_to_cpool(p) + psnshade_to_cpool(p)

      ! get the time step total maintenance respiration
      ! These fluxes should already be in gC/m2/s

      mr = leaf_mr(p) + froot_mr(p)
      if ( iswood(p) ) then
        mr = mr + livestem_mr(p) + livecroot_mr(p)
      else if ( ivt(p) >= npcropmin ) then
        if ( croplive(p) ) mr = mr + livestem_mr(p) + grain_mr(p)
      end if

      ! carbon flux available for allocation
      availc(p) = gpp(p) - mr

      ! new code added for isotope calculations, 7/1/05, PET
      ! If mr > gpp, then some mr comes from gpp, the rest comes from
      ! cpool (xsmr)
      if ( mr > 0._rk8 .and. availc(p) < 0._rk8 ) then
        curmr = gpp(p)
        curmr_ratio = curmr / mr
      else
        curmr_ratio = 1._rk8
      end if
      leaf_curmr(p) = leaf_mr(p) * curmr_ratio
      leaf_xsmr(p) = leaf_mr(p) - leaf_curmr(p)
      froot_curmr(p) = froot_mr(p) * curmr_ratio
      froot_xsmr(p) = froot_mr(p) - froot_curmr(p)
      livestem_curmr(p) = livestem_mr(p) * curmr_ratio
      livestem_xsmr(p) = livestem_mr(p) - livestem_curmr(p)
      livecroot_curmr(p) = livecroot_mr(p) * curmr_ratio
      livecroot_xsmr(p) = livecroot_mr(p) - livecroot_curmr(p)
      grain_curmr(p) = grain_mr(p) * curmr_ratio
      grain_xsmr(p) = grain_mr(p) - grain_curmr(p)

      ! no allocation when available c is negative
      availc(p) = max(availc(p),0.0_rk8)

      ! test for an xsmrpool deficit
      if ( xsmrpool(p) < 0.0_rk8 ) then
        ! Running a deficit in the xsmrpool, so the first priority is to let
        ! some availc from this timestep accumulate in xsmrpool.
        ! Determine rate of recovery for xsmrpool deficit

        xsmrpool_recover(p) = -xsmrpool(p)/(dayscrecover*secspday)
        if ( xsmrpool_recover(p) < availc(p) ) then
          ! available carbon reduced by amount for xsmrpool recovery
          availc(p) = availc(p) - xsmrpool_recover(p)
        else
          ! all of the available carbon goes to xsmrpool recovery
          xsmrpool_recover(p) = availc(p)
          availc(p) = 0.0_rk8
        end if
        cpool_to_xsmrpool(p) = xsmrpool_recover(p)
      end if

      f1 = froot_leaf(ivt(p))
      f2 = croot_stem(ivt(p))

      ! modified wood allocation to be 2.2 at npp=800 gC/m2/yr, 0.2 at npp=0,
      ! constrained so that it does not go lower than 0.2 (under negative
      ! annsum_npp)
      ! This variable allocation is only for trees. Shrubs have a constant
      ! allocation as specified in the pft-physiology file.
      !  The value is also used as a trigger here:
      !        -1.0 means to use the dynamic allocation (trees).
      if (stem_leaf(ivt(p)) == -1._rk8) then
        f3 = (2.7_rk8/(1.0_rk8+exp(-0.004_rk8 * &
               (annsum_npp(p) - 300.0_rk8)))) - 0.4_rk8
      else
        f3 = stem_leaf(ivt(p))
      end if

      f4 = flivewd(ivt(p))
      g1 = grperc(ivt(p))
      g2 = grpnow(ivt(p))
      cnl = leafcn(ivt(p))
      cnfr = frootcn(ivt(p))
      cnlw = livewdcn(ivt(p))
      cndw = deadwdcn(ivt(p))

      ! calculate f1 to f5 for prog crops following AgroIBIS subr phenocrop

      f5 = 0._rk8 ! continued intializations from above

      if ( ivt(p) >= npcropmin ) then ! skip 2 generic crops

        if ( croplive(p) ) then
          ! same phases appear in subroutine CropPhenology

          ! Phase 1 completed:
          ! ==================
          ! if hui is less than the number of gdd needed for filling of grain
          ! leaf emergence also has to have taken place for lai changes to occur
          ! and carbon assimilation
          ! Next phase: leaf emergence to start of leaf decline

          if (leafout(p) >= huileaf(p) .and. hui(p) < huigrain(p)) then

            ! allocation rules for crops based on maturity and linear decrease
            ! of amount allocated to roots over course of the growing season

            if ( peaklai(p) == 1 ) then ! lai at maximum allowed
              arepr(p) = 0._rk8
              aleaf(p) = 1.e-5_rk8
              astem(p) = 0._rk8
              aroot(p) = 1._rk8 - arepr(p) - aleaf(p) - astem(p)
            else
              arepr(p) = 0._rk8
              aroot(p) = max(0._rk8, min(1._rk8, arooti(ivt(p)) -   &
                          (arooti(ivt(p)) - arootf(ivt(p))) *  &
                           min(1._rk8, hui(p)/gddmaturity(p))))
              ! fraction alloc to leaf (from J Norman alloc curve)
              fleaf = fleafi(ivt(p)) * (exp(-bfact(ivt(p))) -         &
                      exp(-bfact(ivt(p))*hui(p)/huigrain(p))) / &
                      (exp(-bfact(ivt(p)))-1.0_rk8)
              aleaf(p) = max(1.e-5_rk8, (1._rk8 - aroot(p)) * fleaf)
              astem(p) = 1._rk8 - arepr(p) - aleaf(p) - aroot(p)
            end if

            ! AgroIBIS included here an immediate adjustment to aleaf & astem
            ! if the predicted lai from the above allocation coefficients
            ! exceeded laimx.
            ! We have decided to live with lais slightly higher than laimx by
            ! enforcing the cap in the following tstep through the peaklai
            ! logic above.

            astemi(p) = astem(p) ! save for use by equations after shift
            aleafi(p) = aleaf(p) ! to reproductive phenology stage begins
            grain_flag(p) = 0._rk8 ! setting to 0 while in phase 2

            ! Phase 2 completed:
            ! ==================
            ! shift allocation either when enough gdd are accumulated or
            ! maximum number of days has elapsed since planting

          else if ( hui(p) >= huigrain(p) ) then

            aroot(p) = max(0._rk8, min(1._rk8, arooti(ivt(p)) - &
                      (arooti(ivt(p)) - arootf(ivt(p))) * &
                      min(1._rk8, hui(p)/gddmaturity(p))))
            if (astemi(p) > astemf(ivt(p))) then
              astem(p) = max(0._rk8, max(astemf(ivt(p)), astem(p) * &
                         (1._rk8 - min((hui(p)-                 &
                         huigrain(p))/((gddmaturity(p)*declfact(ivt(p)))- &
                         huigrain(p)),1._rk8)**allconss(ivt(p)) )))
            end if
            if ( aleafi(p) > aleaff(ivt(p)) ) then
              aleaf(p) = max(1.e-5_rk8, max(aleaff(ivt(p)), aleaf(p) * &
                            (1._rk8 - min((hui(p)-                    &
                            huigrain(p))/((gddmaturity(p)*declfact(ivt(p)))- &
                            huigrain(p)),1._rk8)**allconsl(ivt(p)) )))
            end if

            ! Beth's retranslocation of leafn, stemn, rootn to organ
            ! Filter excess plant N to retransn pool for organ N
            ! Only do one time then hold grain_flag till onset next season

            ! slevis: Will astem ever = astemf exactly?
            ! Beth's response: ...looks like astem can equal astemf under the
            ! right circumstances.
            ! It might be worth a rewrite to capture what I was trying to do,
            ! but the retranslocation for corn and wheat begins at the beginning
            ! of the grain fill stage, but for soybean I was holding it until
            ! after the leaf and stem decline were complete. Looking at
            ! how astem is calculated, once the stem decline is near complete,
            ! astem should (usually) be set to astemf.
            ! The reason for holding off on soybean is that the retranslocation
            ! scheme begins at the beginning of the grain phase, when the leaf
            ! and stem are still growing, but declining.
            ! Since carbon is still getting allocated and now
            ! there is more nitrogen available, the nitrogen can be divertedi
            ! from grain. For corn and wheat the impact was probably enough to
            ! boost productivity, but for soybean the nitrogen was better off
            ! fulfilling the grain fill. It seems that if the peak lai is
            ! reached for soybean though that this would be bypassed
            ! altogether, not the intended outcome. I checked several of my
            ! output files and they all seemed to be going through the
            ! retranslocation loop for soybean - good news.

            if ( ivt(p) /= nsoybean .or. astem(p) == astemf(ivt(p)) ) then
              if ( grain_flag(p) == 0._rk8 ) then
                t1 = 1 / dt
                leafn_to_retransn(p) = t1 * ((leafc(p) / leafcn(ivt(p))) - &
                                             (leafc(p) / fleafcn(ivt(p))))
                livestemn_to_retransn(p) = t1 * ((livestemc(p) / &
                        livewdcn(ivt(p))) - (livestemc(p) / fstemcn(ivt(p))))
                frootn_to_retransn(p) = 0._rk8
                if ( ffrootcn(ivt(p)) > 0._rk8 ) then
                  frootn_to_retransn(p) = t1 * &
                          ((frootc(p) / frootcn(ivt(p))) - &
                          (frootc(p) / ffrootcn(ivt(p))))
                end if
                grain_flag(p) = 1._rk8
              end if
            end if
            arepr(p) = 1._rk8 - aroot(p) - astem(p) - aleaf(p)
          else                   ! pre emergence
            aleaf(p) = 1.e-5_rk8 ! allocation coefficients should be irrelevant
            astem(p) = 0._rk8    ! because crops have no live carbon pools;
            aroot(p) = 0._rk8    ! this applies to this "else" and to the "else"
            arepr(p) = 0._rk8    ! a few lines down
          end if
          f1 = aroot(p) / aleaf(p)
          f3 = astem(p) / aleaf(p)
          f5 = arepr(p) / aleaf(p)
          g1 = 0.25_rk8
        else   ! .not croplive
          f1 = 0._rk8
          f3 = 0._rk8
          f5 = 0._rk8
          g1 = 0.25_rk8
        end if
      end if

      ! based on available C, use constant allometric relationships to
      ! determine N requirements
      if ( iswood(p) ) then
        c_allometry(p) = (1._rk8+g1)*(1._rk8+f1+f3*(1._rk8+f2))
        n_allometry(p) = 1._rk8/cnl + f1/cnfr + (f3*f4*(1._rk8+f2))/cnlw + &
                       (f3*(1._rk8-f4)*(1._rk8+f2))/cndw
      else if (ivt(p) >= npcropmin) then ! skip generic crops
        cng = graincn(ivt(p))
        c_allometry(p) = (1._rk8+g1)*(1._rk8+f1+f5+f3*(1._rk8+f2))
        n_allometry(p) = 1._rk8/cnl + f1/cnfr + f5/cng + &
                (f3*f4*(1._rk8+f2))/cnlw + &
                (f3*(1._rk8-f4)*(1._rk8+f2))/cndw
      else
        c_allometry(p) = 1._rk8+g1+f1+f1*g1
        n_allometry(p) = 1._rk8/cnl + f1/cnfr
      end if
      plant_ndemand(p) = availc(p)*(n_allometry(p)/c_allometry(p))

      ! retranslocated N deployment depends on seasonal cycle of potential GPP
      ! (requires one year run to accumulate demand)

      tempsum_potential_gpp(p) = tempsum_potential_gpp(p) + gpp(p)

      ! Adding the following line to carry max retransn info to CN Annual Update
      tempmax_retransn(p) = max(tempmax_retransn(p),retransn(p))

      ! Beth's code: crops pull from retransn pool only during grain fill;
      !              retransn pool has N from leaves, stems, and roots for
      !              retranslocation

      if ( ivt(p) >= npcropmin .and. grain_flag(p) == 1._rk8 ) then
        avail_retransn(p) = plant_ndemand(p)
      else if ( ivt(p) < npcropmin .and. annsum_potential_gpp(p) > 0._rk8 ) then
        avail_retransn(p) = (annmax_retransn(p)/2._rk8) * &
                (gpp(p)/annsum_potential_gpp(p))/dt
      else
        avail_retransn(p) = 0.0_rk8
      end if

      ! make sure available retrans N doesn't exceed storage
      avail_retransn(p) = min(avail_retransn(p), retransn(p)/dt)

      ! modify plant N demand according to the availability of
      ! retranslocated N
      ! take from retransn pool at most the flux required to meet
      ! plant ndemand

      if (plant_ndemand(p) > avail_retransn(p)) then
        retransn_to_npool(p) = avail_retransn(p)
      else
        retransn_to_npool(p) = plant_ndemand(p)
      end if
      plant_ndemand(p) = plant_ndemand(p) - retransn_to_npool(p)

    end do ! end pft loop

    ! now use the p2c routine to get the column-averaged plant_ndemand
    call p2c(num_soilc,filter_soilc,plant_ndemand,col_plant_ndemand)

#ifndef NITRIF_DENITRIF
    ! column loops to resolve plant/heterotroph competition for mineral N

    ! init sminn_tot
    do fc = 1, num_soilc
      c = filter_soilc(fc)
      sminn_tot(c) = 0.
    end do

    do j = 1, nlevdecomp
      do fc = 1, num_soilc
         c = filter_soilc(fc)
         sminn_tot(c) = sminn_tot(c) + sminn_vr(c,j) * dzsoi_decomp(j)
      end do
    end do

    do j = 1, nlevdecomp
      do fc = 1, num_soilc
        c = filter_soilc(fc)
        if (sminn_tot(c) > 0.) then
          nuptake_prof(c,j) = sminn_vr(c,j) / sminn_tot(c)
        else
          nuptake_prof(c,j) = nfixation_prof(c,j)
        end if
        sum_ndemand_vr(c,j) = col_plant_ndemand(c) * nuptake_prof(c,j) + &
                potential_immob_vr(c,j)
      end do
    end do

    do j = 1, nlevdecomp
      do fc = 1, num_soilc
        c = filter_soilc(fc)
        l = clandunit(c)
        if ( sum_ndemand_vr(c,j)*dt < sminn_vr(c,j) ) then
           ! N availability is not limiting immobilization or plant
           ! uptake, and both can proceed at their potential rates
           nlimit(c,j) = 0
           fpi_vr(c,j) = 1.0_rk8
           actual_immob_vr(c,j) = potential_immob_vr(c,j)
           sminn_to_plant_vr(c,j) = col_plant_ndemand(c) * nuptake_prof(c,j)
        else if ( carbon_only ) then !.or. &
          !       (crop_supln .and. (itypelun(l) == istcrop) .and. &
          !       (ivt(pfti(c)) >= npcropmin)) )then
          ! this code block controls the addition of N to sminn pool
          ! to eliminate any N limitation, when Carbon_Only is set.
          !  This lets the model behave essentially as a carbon-only model,
          ! but with the benefit of keeping track of the N additions needed to
          ! eliminate N limitations, so there is still a diagnostic quantity
          ! that describes the degree of N limitation at steady-state.
          nlimit(c,j) = 1
          fpi_vr(c,j) = 1.0_rk8
          actual_immob_vr(c,j) = potential_immob_vr(c,j)
          sminn_to_plant_vr(c,j) =  col_plant_ndemand(c) * nuptake_prof(c,j)
          supplement_to_sminn_vr(c,j) = sum_ndemand_vr(c,j) - (sminn_vr(c,j)/dt)
        else
          ! N availability can not satisfy the sum of immobilization and
          ! plant growth demands, so these two demands compete for available
          ! soil mineral N resource.

          nlimit(c,j) = 1
          if ( sum_ndemand_vr(c,j) > 0.0_rk8 ) then
            actual_immob_vr(c,j) = (sminn_vr(c,j)/dt) * &
                    (potential_immob_vr(c,j) / sum_ndemand_vr(c,j))
          else
            actual_immob_vr(c,j) = 0.0_rk8
          end if

          if ( potential_immob_vr(c,j) > 0.0_rk8 ) then
            fpi_vr(c,j) = actual_immob_vr(c,j) / potential_immob_vr(c,j)
          else
            fpi_vr(c,j) = 0.0_rk8
          end if
          sminn_to_plant_vr(c,j) = (sminn_vr(c,j)/dt) - actual_immob_vr(c,j)
        end if
      end do
    end do

    ! sum up N fluxes to plant
    do j = 1, nlevdecomp
      do fc = 1, num_soilc
        c = filter_soilc(fc)
        sminn_to_plant(c) = sminn_to_plant(c) + &
                sminn_to_plant_vr(c,j) * dzsoi_decomp(j)
      end do
    end do

    ! give plants a second pass to see if there is any mineral N left over
    ! with which to satisfy residual N demand.
    do fc = 1, num_soilc
      c = filter_soilc(fc)
      residual_sminn(c) = 0._rk8
    end do

    ! sum up total N left over after initial plant and immobilization fluxes
    do fc = 1, num_soilc
      c = filter_soilc(fc)
      residual_plant_ndemand(c) = col_plant_ndemand(c) - sminn_to_plant(c)
    end do
    do j = 1, nlevdecomp
      do fc = 1, num_soilc
        c = filter_soilc(fc)
        if ( residual_plant_ndemand(c) > 0._rk8 ) then
          if ( nlimit(c,j) == 0 ) then
            residual_sminn_vr(c,j) = max(sminn_vr(c,j) - &
                    (actual_immob_vr(c,j) + sminn_to_plant_vr(c,j) ) * dt, 0._rk8)
            residual_sminn(c) = residual_sminn(c) + &
                    residual_sminn_vr(c,j) * dzsoi_decomp(j)
          else
            residual_sminn_vr(c,j)  = 0._rk8
          end if
        end if
      end do
    end do

    ! distribute residual N to plants
    do j = 1, nlevdecomp
      do fc = 1, num_soilc
        c = filter_soilc(fc)
        if ( residual_plant_ndemand(c) > 0._rk8 .and. &
             residual_sminn(c) > 0._rk8 .and. nlimit(c,j) == 0 ) then
          sminn_to_plant_vr(c,j) = sminn_to_plant_vr(c,j) + &
                  residual_sminn_vr(c,j) * &
                  min(( residual_plant_ndemand(c) *  dt ) / &
                  residual_sminn(c), 1._rk8) / dt
        end if
      end do
    end do

    ! re-sum up N fluxes to plant
    do fc = 1, num_soilc
      c = filter_soilc(fc)
      sminn_to_plant(c) = 0._rk8
    end do
    do j = 1, nlevdecomp
      do fc = 1, num_soilc
        c = filter_soilc(fc)
        sminn_to_plant(c) = sminn_to_plant(c) + &
                sminn_to_plant_vr(c,j) * dzsoi_decomp(j)
        sum_ndemand_vr(c,j) = potential_immob_vr(c,j) + sminn_to_plant_vr(c,j)
      end do
    end do

    ! under conditions of excess N, some proportion is assumed to
    ! be lost to denitrification, in addition to the constant
    ! proportion lost in the decomposition pathways
    do j = 1, nlevdecomp
      do fc = 1, num_soilc
        c = filter_soilc(fc)
        if ( (sminn_to_plant_vr(c,j) + actual_immob_vr(c,j))*dt < &
                sminn_vr(c,j) ) then
          sminn_to_denit_excess_vr(c,j) = max(bdnr*((sminn_vr(c,j)/dt) - &
                  sum_ndemand_vr(c,j)),0._rk8)
        else
          sminn_to_denit_excess_vr(c,j) = 0._rk8
        end if
      end do
    end do

    ! sum up N fluxes to immobilization
    do j = 1, nlevdecomp
      do fc = 1, num_soilc
        c = filter_soilc(fc)
        actual_immob(c) = actual_immob(c) + &
                actual_immob_vr(c,j) * dzsoi_decomp(j)
        potential_immob(c) = potential_immob(c) + &
                potential_immob_vr(c,j) * dzsoi_decomp(j)
      end do
    end do

    do fc = 1, num_soilc
      c = filter_soilc(fc)
      ! calculate the fraction of potential growth that can be
      ! acheived with the N available to plants
      if (col_plant_ndemand(c) > 0.0_rk8) then
        fpg(c) = sminn_to_plant(c) / col_plant_ndemand(c)
      else
        fpg(c) = 1.0_rk8
      end if

      ! calculate the fraction of immobilization realized
      ! (for diagnostic purposes)
      if ( potential_immob(c) > 0.0_rk8 ) then
        fpi(c) = actual_immob(c) / potential_immob(c)
      else
        fpi(c) = 1.0_rk8
      end if
    end do

#else
    !----------NITRIF_DENITRIF-------------!
    ! column loops to resolve plant/heterotroph/nitrifier/denitrifier
    ! competition for mineral N

    ! init total mineral N pools
    do fc = 1, num_soilc
      c = filter_soilc(fc)
      sminn_tot(c) = 0.
    end do

    ! sum up total mineral N pools
    do j = 1, nlevdecomp
      do fc = 1, num_soilc
        c = filter_soilc(fc)
        sminn_tot(c) = sminn_tot(c) + &
                (smin_no3_vr(c,j) + smin_nh4_vr(c,j)) * dzsoi_decomp(j)
      end do
    end do

    ! define N uptake profile for initial vertical distribution of plant N
    ! uptake, assuming plant seeks N from where it is most abundant
    do j = 1, nlevdecomp
      do fc = 1, num_soilc
        c = filter_soilc(fc)
        if (sminn_tot(c) > 0.) then
          nuptake_prof(c,j) = sminn_vr(c,j) / sminn_tot(c)
        else
          nuptake_prof(c,j) = nfixation_prof(c,j)
        end if
      end do
    end do

    ! main column/vertical loop
    do j = 1, nlevdecomp
      do fc = 1, num_soilc
        c = filter_soilc(fc)
        l = clandunit(c)

        !  first compete for nh4
        sum_nh4_demand(c,j) = col_plant_ndemand(c) * nuptake_prof(c,j) + &
                potential_immob_vr(c,j) + pot_f_nit_vr(c,j)
        sum_nh4_demand_scaled(c,j) = col_plant_ndemand(c)*nuptake_prof(c,j) * &
                compet_plant_nh4 + potential_immob_vr(c,j)*compet_decomp_nh4 + &
                pot_f_nit_vr(c,j)*compet_nit
        if ( sum_nh4_demand(c,j)*dt < smin_nh4_vr(c,j) ) then
          ! NH4 availability is not limiting immobilization or plant
          ! uptake, and all can proceed at their potential rates
          nlimit_nh4(c,j) = 0
          fpi_nh4_vr(c,j) = 1.0_rk8
          actual_immob_nh4_vr(c,j) = potential_immob_vr(c,j)
          smin_nh4_to_plant_vr(c,j) = col_plant_ndemand(c) * nuptake_prof(c,j)

          f_nit_vr(c,j) = pot_f_nit_vr(c,j)

        else
          ! NH4 availability can not satisfy the sum of immobilization,
          ! nitrification, and plant growth demands, so these three demands
          ! compete for available soil mineral NH4 resource.
          nlimit_nh4(c,j) = 1
          if ( sum_nh4_demand(c,j) > 0.0_rk8 ) then
            actual_immob_nh4_vr(c,j) = min((smin_nh4_vr(c,j)/dt)*&
                 (potential_immob_vr(c,j)*compet_decomp_nh4 / &
                 sum_nh4_demand_scaled(c,j)), potential_immob_vr(c,j))
            smin_nh4_to_plant_vr(c,j) = min((smin_nh4_vr(c,j)/dt)*&
                 (col_plant_ndemand(c)*nuptake_prof(c,j)*compet_plant_nh4 / &
                 sum_nh4_demand_scaled(c,j)), &
                 col_plant_ndemand(c)*nuptake_prof(c,j))
            f_nit_vr(c,j) =  min((smin_nh4_vr(c,j)/dt)*(pot_f_nit_vr(c,j)*&
               compet_nit / sum_nh4_demand_scaled(c,j)), pot_f_nit_vr(c,j))
          else
            actual_immob_nh4_vr(c,j) = 0.0_rk8
            smin_nh4_to_plant_vr(c,j) = 0.0_rk8
            f_nit_vr(c,j) = 0.0_rk8
          end if

          if ( potential_immob_vr(c,j) > 0.0_rk8 ) then
            fpi_nh4_vr(c,j) = actual_immob_nh4_vr(c,j) / potential_immob_vr(c,j)
          else
            fpi_nh4_vr(c,j) = 0.0_rk8
          end if
        end if

        ! next compete for no3
        sum_no3_demand(c,j) = (col_plant_ndemand(c)*nuptake_prof(c,j)-&
           smin_nh4_to_plant_vr(c,j)) + (potential_immob_vr(c,j)-&
           actual_immob_nh4_vr(c,j)) + pot_f_denit_vr(c,j)
        sum_no3_demand_scaled(c,j) = (col_plant_ndemand(c)*&
           nuptake_prof(c,j)-smin_nh4_to_plant_vr(c,j))*compet_plant_no3 + &
              (potential_immob_vr(c,j)-actual_immob_nh4_vr(c,j))*&
              compet_decomp_no3 + pot_f_denit_vr(c,j)*compet_denit

        if ( sum_no3_demand(c,j)*dt < smin_no3_vr(c,j) ) then
          ! NO3 availability is not limiting immobilization or plant
          ! uptake, and all can proceed at their potential rates
          nlimit_no3(c,j) = 1
          fpi_no3_vr(c,j) = 1.0_rk8 -  fpi_nh4_vr(c,j)
            actual_immob_no3_vr(c,j) = (potential_immob_vr(c,j) - &
                    actual_immob_nh4_vr(c,j))
          smin_no3_to_plant_vr(c,j) = (col_plant_ndemand(c) * &
                  nuptake_prof(c,j) - smin_nh4_to_plant_vr(c,j))
          f_denit_vr(c,j) = pot_f_denit_vr(c,j)
        else
          ! NO3 availability can not satisfy the sum of immobilization,
          ! denitrification, and plant growth demands, so these three
          ! demands compete for available soil mineral NO3 resource.
          nlimit_no3(c,j) = 1
          if ( sum_no3_demand(c,j) > 0.0_rk8 ) then
            actual_immob_no3_vr(c,j) = min((smin_no3_vr(c,j)/dt)*&
                 ((potential_immob_vr(c,j)-actual_immob_nh4_vr(c,j))*&
                 compet_decomp_no3 / sum_no3_demand_scaled(c,j)), &
                    potential_immob_vr(c,j)-actual_immob_nh4_vr(c,j))
            smin_no3_to_plant_vr(c,j) = min((smin_no3_vr(c,j)/dt)*&
                 ((col_plant_ndemand(c)*nuptake_prof(c,j)-&
                 smin_nh4_to_plant_vr(c,j))*compet_plant_no3 / &
                 sum_no3_demand_scaled(c,j)), &
                 col_plant_ndemand(c)*nuptake_prof(c,j) - &
                 smin_nh4_to_plant_vr(c,j))
            f_denit_vr(c,j) =  min((smin_no3_vr(c,j)/dt)*&
                 (pot_f_denit_vr(c,j)*compet_denit / &
                 sum_no3_demand_scaled(c,j)), pot_f_denit_vr(c,j))
          else
            actual_immob_no3_vr(c,j) = 0.0_rk8
            smin_no3_to_plant_vr(c,j) = 0.0_rk8
            f_denit_vr(c,j) = 0.0_rk8
          end if

          if ( potential_immob_vr(c,j) > 0.0_rk8 ) then
            fpi_no3_vr(c,j) = actual_immob_no3_vr(c,j) / potential_immob_vr(c,j)
          else
            fpi_no3_vr(c,j) = 0.0_rk8
          end if

        end if

        ! n2o emissions: n2o from nitr is const fraction, n2o from denitr
        ! is calculated in nitrif_denitrif
        f_n2o_nit_vr(c,j) = f_nit_vr(c,j) * nitrif_n2o_loss_frac
        f_n2o_denit_vr(c,j) = f_denit_vr(c,j) / &
                (1._rk8 + n2_n2o_ratio_denit_vr(c,j))


        ! samy : Calculation of total N2O emission from soil to atmosphere (
        ! retrieved from Dynamical Land Ecosystem Model (DLEM))
        f_n2o_tot_vr(c,j) = 0.021_rk8 * &
          ((0.001_rk8 * f_nit_vr(c,j)) + f_denit_vr(c,j)) * &
          ((10_rk8**(h2osoi_vol(c,j)/fr_WFPS(c,j))*0.026_rk8 - 1.66_rk8) / &
          (1_rk8 + 10_rk8**(h2osoi_vol(c,j)/fr_WFPS(c,j)))*0.026_rk8 - 1.66_rk8)

        ! this code block controls the addition of N to sminn pool
        ! to eliminate any N limitation, when Carbon_Only is set.  This lets the
        ! model behave essentially as a carbon-only model, but with the
        ! benefit of keeping track of the N additions needed to
        ! eliminate N limitations, so there is still a diagnostic quantity
        ! that describes the degree of N limitation at steady-state.

        if ( carbon_only ) then !.or. &
!            (crop_supln .and. (itypelun(l) == istcrop) .and. &
!            (ivt(pfti(c)) >= npcropmin)) ) then

          if ( fpi_no3_vr(c,j) + fpi_nh4_vr(c,j) < 1._rk8 ) then
            fpi_nh4_vr(c,j) = 1.0_rk8 - fpi_no3_vr(c,j)
            supplement_to_sminn_vr(c,j) = (potential_immob_vr(c,j) - &
                    actual_immob_no3_vr(c,j)) - actual_immob_nh4_vr(c,j)
            ! update to new values that satisfy demand
            actual_immob_nh4_vr(c,j) = potential_immob_vr(c,j) - &
                    actual_immob_no3_vr(c,j)
          end if
          if ( smin_no3_to_plant_vr(c,j) + smin_nh4_to_plant_vr(c,j) < &
                  col_plant_ndemand(c)*nuptake_prof(c,j) ) then
            ! use old values
            supplement_to_sminn_vr(c,j) = supplement_to_sminn_vr(c,j) + &
                 (col_plant_ndemand(c)*nuptake_prof(c,j) - &
                 smin_no3_to_plant_vr(c,j)) - smin_nh4_to_plant_vr(c,j)
            smin_nh4_to_plant_vr(c,j) = col_plant_ndemand(c) * &
                    nuptake_prof(c,j) - smin_no3_to_plant_vr(c,j)
          end if
          sminn_to_plant_vr(c,j) = smin_no3_to_plant_vr(c,j) + &
                  smin_nh4_to_plant_vr(c,j)
        end if

        ! sum up no3 and nh4 fluxes
        fpi_vr(c,j) = fpi_no3_vr(c,j) + fpi_nh4_vr(c,j)
        sminn_to_plant_vr(c,j) = smin_no3_to_plant_vr(c,j) + &
                smin_nh4_to_plant_vr(c,j)
        actual_immob_vr(c,j) = actual_immob_no3_vr(c,j) + &
                actual_immob_nh4_vr(c,j)
      end do
    end do

    do fc = 1, num_soilc
      c = filter_soilc(fc)
      ! sum up N fluxes to plant after initial competition
      sminn_to_plant(c) = 0._rk8
    end do
    do j = 1, nlevdecomp
      do fc = 1, num_soilc
        c = filter_soilc(fc)
        sminn_to_plant(c) = sminn_to_plant(c) + &
                sminn_to_plant_vr(c,j) * dzsoi_decomp(j)
      end do
    end do

    ! give plants a second pass to see if there is any mineral N left over
    ! with which to satisfy residual N demand.
    ! first take frm nh4 pool; then take from no3 pool
    do fc = 1, num_soilc
      c = filter_soilc(fc)
      residual_plant_ndemand(c) = col_plant_ndemand(c) - sminn_to_plant(c)
      residual_smin_nh4(c) = 0._rk8
    end do
    do j = 1, nlevdecomp
      do fc = 1, num_soilc
        c = filter_soilc(fc)
        if (residual_plant_ndemand(c) > 0._rk8 ) then
          if (nlimit_nh4(c,j) == 0) then
            residual_smin_nh4_vr(c,j) = max(smin_nh4_vr(c,j) - &
                    (actual_immob_vr(c,j) + smin_nh4_to_plant_vr(c,j) ) * &
                    dt, 0._rk8)
            residual_smin_nh4(c) = residual_smin_nh4(c) + &
                    residual_smin_nh4_vr(c,j) * dzsoi_decomp(j)
          else
            residual_smin_nh4_vr(c,j)  = 0._rk8
          end if
          if ( residual_smin_nh4(c) > 0._rk8 .and. nlimit_nh4(c,j) == 0 ) then
            smin_nh4_to_plant_vr(c,j) = smin_nh4_to_plant_vr(c,j) + &
                    residual_smin_nh4_vr(c,j) * &
                    min(( residual_plant_ndemand(c) * dt ) / &
                    residual_smin_nh4(c), 1._rk8) / dt
          end if
        end if
      end do
    end do

    ! re-sum up N fluxes to plant after second pass for nh4
    do fc = 1, num_soilc
      c = filter_soilc(fc)
      sminn_to_plant(c) = 0._rk8
    end do
    do j = 1, nlevdecomp
      do fc = 1, num_soilc
        c = filter_soilc(fc)
        sminn_to_plant_vr(c,j) = smin_nh4_to_plant_vr(c,j) + &
                smin_no3_to_plant_vr(c,j)
        sminn_to_plant(c) = sminn_to_plant(c) + &
                (sminn_to_plant_vr(c,j)) * dzsoi_decomp(j)
      end do
    end do

    !
    ! and now do second pass for no3
    do fc=1,num_soilc
      c = filter_soilc(fc)
      residual_plant_ndemand(c) = col_plant_ndemand(c) - sminn_to_plant(c)
      residual_smin_no3(c) = 0._rk8
    end do

    do j = 1, nlevdecomp
      do fc = 1, num_soilc
        c = filter_soilc(fc)
        if (residual_plant_ndemand(c) > 0._rk8 ) then
          if ( nlimit_no3(c,j) == 0 ) then
            residual_smin_no3_vr(c,j) = max(smin_no3_vr(c,j) - &
                    (actual_immob_vr(c,j) + smin_no3_to_plant_vr(c,j) ) * &
                    dt, 0._rk8)
            residual_smin_no3(c) = residual_smin_no3(c) + &
                    residual_smin_no3_vr(c,j) * dzsoi_decomp(j)
          else
            residual_smin_no3_vr(c,j)  = 0._rk8
          end if

          if ( residual_smin_no3(c) > 0._rk8 .and. nlimit_no3(c,j) == 0 ) then
            smin_no3_to_plant_vr(c,j) = smin_no3_to_plant_vr(c,j) + &
                    residual_smin_no3_vr(c,j) * &
                    min(( residual_plant_ndemand(c) *  dt ) / &
                    residual_smin_no3(c), 1._rk8) / dt
          end if
        end if
      end do
    end do

    ! re-sum up N fluxes to plant after second passes of both no3 and nh4
    do fc=1,num_soilc
      c = filter_soilc(fc)
      sminn_to_plant(c) = 0._rk8
    end do
    do j = 1, nlevdecomp
      do fc = 1, num_soilc
        c = filter_soilc(fc)
        sminn_to_plant_vr(c,j) = smin_nh4_to_plant_vr(c,j) + &
                smin_no3_to_plant_vr(c,j)
        sminn_to_plant(c) = sminn_to_plant(c) + &
                (sminn_to_plant_vr(c,j)) * dzsoi_decomp(j)
      end do
    end do

    ! sum up N fluxes to immobilization
    do fc = 1, num_soilc
      c = filter_soilc(fc)
      actual_immob(c) = 0._rk8
      potential_immob(c) = 0._rk8
    end do
    do j = 1, nlevdecomp
      do fc = 1, num_soilc
        c = filter_soilc(fc)
        actual_immob(c) = actual_immob(c) + &
                actual_immob_vr(c,j) * dzsoi_decomp(j)
        potential_immob(c) = potential_immob(c) + &
                potential_immob_vr(c,j) * dzsoi_decomp(j)
      end do
    end do


    do fc = 1, num_soilc
      c = filter_soilc(fc)
      ! calculate the fraction of potential growth that can be
      ! acheived with the N available to plants
      if ( col_plant_ndemand(c) > 0.0_rk8 ) then
        fpg(c) = sminn_to_plant(c) / col_plant_ndemand(c)
      else
        fpg(c) = 1._rk8
      end if

      ! calculate the fraction of immobilization realized
      ! (for diagnostic purposes)
      if ( potential_immob(c) > 0.0_rk8 ) then
        fpi(c) = actual_immob(c) / potential_immob(c)
      else
        fpi(c) = 1._rk8
      end if
    end do ! end of column loops
#endif

    ! start new pft loop to distribute the available N between the
    ! competing pfts on the basis of relative demand, and allocate C and N to
    ! new growth and storage

    do fp = 1, num_soilp
      p = filter_soilp(fp)
      c = pcolumn(p)

      ! set some local allocation variables
      f1 = froot_leaf(ivt(p))
      f2 = croot_stem(ivt(p))

      ! modified wood allocation to be 2.2 at npp=800 gC/m2/yr, 0.2 at npp=0,
      ! constrained so that it does not go lower than 0.2 (under negative
      ! annsum_npp)
      ! There was an error in this formula in previous version, where the
      ! coefficient was 0.004 instead of 0.0025.
      ! This variable allocation is only for trees. Shrubs have a constant
      ! allocation as specified in the pft-physiology file.
      !  The value is also used as a trigger here:
      !         -1.0 means to use the dynamic allocation (trees).
      if ( stem_leaf(ivt(p)) == -1._rk8 ) then
        f3 = (2.7_rk8/(1.0_rk8 + &
          exp(-0.004_rk8*(annsum_npp(p) - 300.0_rk8)))) - 0.4_rk8
      else
        f3 = stem_leaf(ivt(p))
      end if

      f4 = flivewd(ivt(p))
      g1 = grperc(ivt(p))
      g2 = grpnow(ivt(p))
      cnl = leafcn(ivt(p))
      cnfr = frootcn(ivt(p))
      cnlw = livewdcn(ivt(p))
      cndw = deadwdcn(ivt(p))
      fcur = fcur2(ivt(p))

      if ( ivt(p) >= npcropmin ) then ! skip 2 generic crops
        if ( croplive(p) ) then
          f1 = aroot(p) / aleaf(p)
          f3 = astem(p) / aleaf(p)
          f5 = arepr(p) / aleaf(p)
          g1 = 0.25_rk8
        else
          f1 = 0._rk8
          f3 = 0._rk8
          f5 = 0._rk8
          g1 = 0.25_rk8
        end if
      end if

      ! increase fcur linearly with ndays_active, until fcur reaches 1.0 at
      ! ndays_active = days/year.
      ! This prevents the continued storage of C and N.
      ! turning off this correction (PET, 12/11/03), instead using bgtr in
      ! phenology algorithm.
      !fcur = fcur + (1._rk8 - fcur)*lgsf(p)
      sminn_to_npool(p) = plant_ndemand(p) * fpg(c)
      plant_nalloc(p) = sminn_to_npool(p) + retransn_to_npool(p)

      ! calculate the associated carbon allocation, and the excess
      ! carbon flux that must be accounted for through downregulation
      plant_calloc(p) = plant_nalloc(p) * (c_allometry(p)/n_allometry(p))

      ! samy : the following condition applied to obtain reasonable values
      ! of downregulation

      !if (ivt(p) == nbrdlf_evr_trp_tree) then
      !  plant_calloc(p) = 0.9_rk8 * plant_nalloc(p) * &
      !                (c_allometry(p)/n_allometry(p))
      !else
      !  plant_calloc(p) = 0.8_rk8 * plant_nalloc(p) * &
      !                (c_allometry(p)/n_allometry(p))
      !end if

      excess_cflux(p) = availc(p) - plant_calloc(p)

      ! reduce gpp fluxes due to N limitation
      if ( gpp(p) > 0.0_rk8 ) then
        downreg(p) = excess_cflux(p)/gpp(p)
        psnsun_to_cpool(p) = psnsun_to_cpool(p)*(1._rk8 - downreg(p))
        psnshade_to_cpool(p) = psnshade_to_cpool(p)*(1._rk8 - downreg(p))
        if ( use_c13 ) then
          c13_psnsun_to_cpool(p) = c13_psnsun_to_cpool(p)*(1._rk8 - downreg(p))
          c13_psnshade_to_cpool(p) = &
                  c13_psnshade_to_cpool(p)*(1._rk8 - downreg(p))
        end if
        if ( use_c14 ) then
          c14_psnsun_to_cpool(p) = c14_psnsun_to_cpool(p)*(1._rk8 - downreg(p))
          c14_psnshade_to_cpool(p) = &
                  c14_psnshade_to_cpool(p)*(1._rk8 - downreg(p))
        end if
      end if

      ! calculate the amount of new leaf C dictated by these allocation
      ! decisions, and calculate the daily fluxes of C and N to current
      ! growth and storage pools

      ! fcur is the proportion of this day's growth that is displayed now,
      ! the remainder going into storage for display next year through the
      ! transfer pools

      nlc = plant_calloc(p) / c_allometry(p)
      cpool_to_leafc(p) = nlc * fcur
      cpool_to_leafc_storage(p) = nlc * (1._rk8 - fcur)
      cpool_to_frootc(p) = nlc * f1 * fcur
      cpool_to_frootc_storage(p) = nlc * f1 * (1._rk8 - fcur)
      if ( iswood(p) ) then
        cpool_to_livestemc(p) = nlc * f3 * f4 * fcur
        cpool_to_livestemc_storage(p) = nlc * f3 * f4 * (1._rk8 - fcur)
        cpool_to_deadstemc(p) = nlc * f3 * (1._rk8 - f4) * fcur
        cpool_to_deadstemc_storage(p)  = nlc * f3 * (1._rk8 - f4) * (1._rk8 - fcur)
        cpool_to_livecrootc(p) = nlc * f2 * f3 * f4 * fcur
        cpool_to_livecrootc_storage(p) = nlc * f2 * f3 * f4 * (1._rk8 - fcur)
        cpool_to_deadcrootc(p) = nlc * f2 * f3 * (1._rk8 - f4) * fcur
        cpool_to_deadcrootc_storage(p) = nlc * f2 * f3 * (1._rk8 - f4) * &
                                         (1._rk8 - fcur)
      end if
      if ( ivt(p) >= npcropmin ) then ! skip 2 generic crops
        cpool_to_livestemc(p) = nlc * f3 * f4 * fcur
        cpool_to_livestemc_storage(p) = nlc * f3 * f4 * (1._rk8 - fcur)
        cpool_to_deadstemc(p) = nlc * f3 * (1._rk8 - f4) * fcur
        cpool_to_deadstemc_storage(p) = nlc * f3 * (1._rk8 - f4) * (1._rk8 - fcur)
        cpool_to_livecrootc(p) = nlc * f2 * f3 * f4 * fcur
        cpool_to_livecrootc_storage(p) = nlc * f2 * f3 * f4 * (1._rk8 - fcur)
        cpool_to_deadcrootc(p) = nlc * f2 * f3 * (1._rk8 - f4) * fcur
        cpool_to_deadcrootc_storage(p) = nlc * f2 * f3 * (1._rk8 - f4) * &
                                         (1._rk8 - fcur)
        cpool_to_grainc(p) = nlc * f5 * fcur
        cpool_to_grainc_storage(p) = nlc * f5 * (1._rk8 -fcur)
      end if

      ! corresponding N fluxes
      npool_to_leafn(p) = (nlc / cnl) * fcur
      npool_to_leafn_storage(p) = (nlc / cnl) * (1._rk8 - fcur)
      npool_to_frootn(p) = (nlc * f1 / cnfr) * fcur
      npool_to_frootn_storage(p) = (nlc * f1 / cnfr) * (1._rk8 - fcur)
      if ( iswood(p) ) then
        npool_to_livestemn(p) = (nlc * f3 * f4 / cnlw) * fcur
        npool_to_livestemn_storage(p)  = (nlc * f3 * f4 / cnlw) * (1._rk8 - fcur)
        npool_to_deadstemn(p) = (nlc * f3 * (1._rk8 - f4) / cndw) * fcur
        npool_to_deadstemn_storage(p)  = (nlc * f3 * (1._rk8 - f4) / cndw) * &
                                         (1._rk8 - fcur)
        npool_to_livecrootn(p) = (nlc * f2 * f3 * f4 / cnlw) * fcur
        npool_to_livecrootn_storage(p) = (nlc * f2 * f3 * f4 / cnlw) * &
                                         (1._rk8 - fcur)
        npool_to_deadcrootn(p) = (nlc * f2 * f3 * (1._rk8 - f4) / cndw) * fcur
        npool_to_deadcrootn_storage(p) = (nlc * f2 * f3 * &
                (1._rk8 - f4) / cndw) * (1._rk8 - fcur)
      end if
      if ( ivt(p) >= npcropmin ) then ! skip 2 generic crops
        cng = graincn(ivt(p))
        npool_to_livestemn(p) = (nlc * f3 * f4 / cnlw) * fcur
        npool_to_livestemn_storage(p) = (nlc * f3 * f4 / cnlw) * (1._rk8 - fcur)
        npool_to_deadstemn(p) = (nlc * f3 * (1._rk8 - f4) / cndw) * fcur
        npool_to_deadstemn_storage(p) = (nlc * f3 * (1._rk8 - f4) / cndw) * &
                                        (1._rk8 - fcur)
        npool_to_livecrootn(p) = (nlc * f2 * f3 * f4 / cnlw) * fcur
        npool_to_livecrootn_storage(p) = (nlc * f2 * f3 * f4 / cnlw) * &
                                         (1._rk8 - fcur)
        npool_to_deadcrootn(p) = (nlc * f2 * f3 * (1._rk8 - f4) / cndw) * fcur
        npool_to_deadcrootn_storage(p) = (nlc * f2 * f3 * &
                (1._rk8 - f4) / cndw) * (1._rk8 - fcur)
        npool_to_grainn(p) = (nlc * f5 / cng) * fcur
        npool_to_grainn_storage(p) = (nlc * f5 / cng) * (1._rk8 -fcur)
      end if

      ! Calculate the amount of carbon that needs to go into growth
      ! respiration storage to satisfy all of the storage growth demands.
      ! Allows for the fraction of growth respiration that is released at the
      ! time of fixation, versus the remaining fraction that is stored for
      ! release at the time of display. Note that all the growth respiration
      ! fluxes that get released on a given timestep are calculated in
      ! growth_resp(), but that the storage of C for growth resp during
      ! display of transferred growth is assigned here.

      gresp_storage = cpool_to_leafc_storage(p) + cpool_to_frootc_storage(p)
      if ( iswood(p) ) then
        gresp_storage = gresp_storage + cpool_to_livestemc_storage(p)
        gresp_storage = gresp_storage + cpool_to_deadstemc_storage(p)
        gresp_storage = gresp_storage + cpool_to_livecrootc_storage(p)
        gresp_storage = gresp_storage + cpool_to_deadcrootc_storage(p)
      end if
      if (ivt(p) >= npcropmin) then ! skip 2 generic crops
        gresp_storage = gresp_storage + cpool_to_livestemc_storage(p)
        gresp_storage = gresp_storage + cpool_to_grainc_storage(p)
      end if
      cpool_to_gresp_storage(p) = gresp_storage * g1 * (1._rk8 - g2)

    end do ! end pft loop

  end subroutine CNAllocation

#endif

end module mod_clm_cnallocation
! vim: tabstop=8 expandtab shiftwidth=2 softtabstop=2
